EPJ manuscript No. 

(will be inserted by the editor) 



o 

IX, 



o 



00 

in 



o 



Cell motility: a viscous fingering analysis of active gels 

M. Ben Amar 1 , 0. V. Manyuhina 1 , and G. Napoli 2 

1 Laboratoire de Physique Statistique, Ecole Normale Superieure, UPMC Univ Paris 06, Universite Paris Diderot, CNRS, 24 rue 
Lhomond, 75005 Paris, France 

2 Dipartimento di Ingegneria dell'Innovazione Universita del Salento, Via per Monteroni - Edificio "Corpo O", 1-73100 Lecce, 
Italy 

Received: date / Revised version: date 



Abstract. The symmetry breaking of the actin network from radial to longitudinal symmetry has been 
identified as the major mechanism for keratocytes (fish cells) motility on solid substrate. For strong friction 
oo ■ coefficient, the two dimensional actin flow which includes the polymerisation at the edge and depolymeri- 

sation in the bulk can be modelled as a Darcy flow, the cell shape and dynamics being then modelled by 
standard complex analysis methods. We use the theory of active gels to describe the orientational order 
of the filaments which varies from the border to the bulk. We show analytically that the reorganisation of 
the cortex is enough to explain the motility of the cell and find the velocity as a function of the orientation 
order parameter in the bulk. 

c3 ' 

PACS. 87.I7.Jj Cell locomotion; chemotaxis - 87.17.Rt Cell adhesion and cell mechanics - 47.20.Gv Vis- 
I 1 cous and viscoelastic instabilities - 47.20.Ky Nonlinearity, bifurcation, and symmetry breaking - 47. 20. Ma 

"""O . Interfacial instabilities 
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O l 1 Introduction 

Three polymers, actin filaments, microtubules and intermediate filaments, are responsible of the eukaryotic cell in- 
tegrity [Tj. Together with molecular motors, these filaments accomplish tasks necessary to the cell survival and cell 
mitosis. Microtubules are known to transport intracellular molecules, chromosomes and organelles inside the cells. They 
play also an important role in multi-cellular organisms during embryogenesis. As an example, for C-elegans at some 
stage of the development, they modify the shape of the embryo from an ovoid to a worm-like shape by redistributing 
the stress originating from the contraction of the actin-myosin network inside the epithelial cells [2J. Assembly and 
disassembly of these filaments [3] (microtubule and actin filaments, see fig. [T]A.) produce forces, generate tension and 
are responsible of the cell shape deformations. 

Here we focus on the crawling of cells like keratocytes (fish epithelial cell) or fibroblasts. By occurring on a solid 
substrate, in vitro, this process induces the protrusion of a lamellipodium mainly due to the polymerisation of the actin 
network in the direction of the displacement. Observations in vitro reveal four well identified stages: the formation of 
this lamellipodium, the adhesion to the substratum, the retraction of the rear, then the de-adhesion. They also clearly 
indicate a shape bifurcation of the whole fiat cellular shape from a static position where the cell is mostly round to a 
C$ ■ crescent shape which accompanies the displacement. Visualization of the actin cortex via fluorescent speckle microscopy 
shows a reorganization of the actin network inside the cell in concordance with this shape symmetry breaking. For 
keratocytes, F-actin moves from the leading edge to the cell core at a rate about 25-60 nm/s. The flux is known to be 
faster at the periphery, decreasing in the center. In motile keratocytes (see fig. QJ3), F-actin moves rewards slowly as 
the cell moves rapidly forwards. In this case the F-actin flow decreases to 10 nm/s. The biological processes involved in 
cell motility are complex and many actors play a decisive role in this migration process like molecular motors (mostly 
myosin) responsible of the deformation of the actin cortex or the Wiskott-Aldrich proteins (WASP) which regulate 
the polymerization process [I] at the level of the lipid membrane. 

Our goal is not to enter into all the details of the microscopic biological mechanisms but to explain the existence 
of a displacement when the cortex of actin is reorganized by the molecular motors. We adopt a macroscopic viewpoint 
by treating the F-actin network as a nematic two-dimensional gel put in a situation out of equilibrium by ATP 
hydrolysis. We will extend the approach of Callan et al. [5] limited to the circular description of the actin network so 
to static cell shapes by taking into account the evolution of the orientation of the filaments as the cell moves. This 
requires a modification of the constitutive equations of active gels, including nonlinearities in the free-energy to explain 
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Fig. 1. A) Mouse fibroblasts stained for actin (blue), microtubules (yellow), and nuclei (green). B) Migrating goldfish skin 
keratocytes stained for actin (blue) and adhesion sites (red). Courtesy Professor Torsten Wittmann, University of California 
http: / /www. ucsf.edu / science-cafe /conversations /wittmann / ) . 

the localization of the cortex at the cell border. The tensorial resulting balance equations are coupled to boundary 
conditions resulting from mechanical equilibrium and observation of the actin network. As an example, it is known 
that the filaments approach the lipid membrane quite tangentially (with some small angle of order twenty degrees [T| 
[5]). The mathematical solutions of these equations with these boundary conditions lead to a free-boundary problem 
where the cell shape cannot be fixed a priori. Treating the disorientated cortex at the cell border as a boundary layer 
in both cases (static and motile) with an isotropic core (static case) or a fully oriented core (motile case), we succeed 
to formulate mathematically the shape equation with the Schwarz function [5J, a technique used for shape drops in 
hydrodynamics and wetting [7]. The actin cortex, once treated as a boundary layer, induces a modification of the 
boundary conditions for the hydrodynamic flow. In particular we explain the destabilization of the static cell for an 
increase of tensile activity of the molecular motors and the motility of the cell, after reorganization of the core for cells 
which remain quasi-circular. 

The paper is organized as follows. First we show that the process of polymerization-depolymerisation induces a flow 
which can be approximated by some Darcy flow when the friction on the substrate is important and we establish the 
boundary conditions of Neumann and Dirichlet type at the origin of the free-boundary problem. In section 3 we present 
the complex analysis formulation [BJ and prove mathematically that the Callan et al [5] approach cannot explain the 
cell motility. Section 4 is a reminder of the theory of active gels [81IWTU] with application to the polar geometry (2D 
cylindrical geometry). Section 5 considers the case where the actin cortex is limited to a boundary layer at the cell 
periphery with a melted and oriented core. Section 6 establishes the modification of the boundary conditions while in 
section 7 we show that at leading order our theory is enough to explain motility for a circular cell with oriented core. 



2 Formulation of the problem 
2.1 The flow origin 

We assume that the shape of the drop is fixed by the actin flow in interaction with the lipid membrane [5 . It is 
the result of the depolymerisation of the actin filaments which is isotropically distributed in the bulk (called fl in 
the following) and of the polymerisation which occurs at the boundary (dfT). Such an assumption is inspired by the 
experimental observations [11] : a flow of actin is visualized and reveals that it is directed towards some center in the 
cell. Such assumption is consistent with the fact that circular cells of constant thickness do not move on the substrate. 
Indeed, we assume that the actin filaments arrive at the cell border with some orientation with the normal and touch 
it with the pole plus. So polymerisation occurs at the border which becomes a sink of actin monomers contrary to the 
bulk which is a source. The conservation of mass is represented by the continuity equation: 

^ + V • Gov) = -pk d + P TV p n5{(r - r mt ) ■ n), (1) 

where p is the actin filament concentration, moving with a local velocity v . The actin gel is assumed incompressible 
and there is a conservation of the total number of actin monomers (either free or organised in filaments or bundles) . p 
is assumed constant inside the cell. Indeed, since our model is two-dimensional and since there exists some thickness 
variation especially at the border this hypothesis may be questionable but we will discard this thickness effect 
and keep a purely two-dimensional model. We call V p the polymerisation velocity, Ti nt the position of a point on the 
interface and n the normal at this point to the interface. A dimensionless shape factor T indicates that the density 
of actin filaments at the border may vary from point to point. It is typically a function of the local curvature or its 
even derivatives. We assume that this shape factor J- is normalized to 1 for the circle which has a constant curvature. 
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We can tune T if we have biological informations. For example Yam et al [llj mention that when cells become motile, 
they are no longer circular and the concentration of actin touching the border is higher at the front than at the rear. 
By taking into account a small volume element centered along the cell border and integrating this equation, we derive 

v n = v • n = C/e y • n - TV P , (2) 

In the bulk (J?), depolymerisation of actin filaments gives 

V-v = -fc rf . (3) 

Moreover, mass conservation of actin filaments imposes that: 

k d / dxdy = V p Tds, (4) 
Jn Jsn 

where s is a curvilinear coordinate along the contour 8Q. As an example, for a circular cell of radius R, we get 
k d R = Vp/2, which coincides with the result derived in [Fj. Note that this mass conservation has to be imposed at each 
step of the cell evolution. It is a constraint which must be satisfied by the shape solution. For non-steady dynamics 
of the motility process, for transitory states or stability studies for example, eq. @ is easily modified by introducing 
the local velocity of the boundary V^^ so we get 

v • n + TV P = V sn ■ n. (5) 



2.2 The flow equation 

We neglect here inertia and the Stokes equation inside the cell is given by 

- VP - £v + ^Av = 0, (6) 

if the cell is assumed to be isotropic and homogeneous. Indeed, the cell contains filaments with molecular motors 
able to transmit stresses which can modify eq. ©. When these filaments are completely disoriented, we will show 
in section 21 that eq. © is not modified although the system is active. It means that inside the cell, there exists 
permanent stresses induced chemically by ATP hydrolysis. When anisotropy of the actin network will be introduced, 
eq. © will be modified (see section 6). In the limit when the friction coefficient £ is strong, we can neglect the viscosity 
r\ and recover from eq. (JBJ) the Darcy law for the flow inside the drop. First we consider the full problem, then the 
limiting cases. 

The flow velocity during motility initiation has been visualised experimentally using multiframe correlation track- 
ing. Initially radial and centripetal, it becomes oriented along the axis of displacement. To satisfy eq. ©, we take 
into account these two cases. We define the stream function ^ which allows to derive the incompressible part of the 
velocity field. In a Cartesian coordinate system, we have 

k d , M k d dV 

for a flow field mostly radial, and for a flow field mainly oriented along the y direction 

v x = w Vy = -—- kd y + Vo . (8) 
being Vq a constant divergence-less velocity. The Stokes equation eq. ([5]) is then transformed into: 

) (9) 
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<k d y + £,V . (10) 

Once equation ([9]) is differentiated with respect to y and eq. (flQ)) with respect to x, the pressure P can be eliminated 
and we get 

r]^ 2 <P - £A<f = 0. (11) 

When the active stress cr° ct = —(8/^6^ is introduced isotropically, it modifies the boundary conditions. The material 
coefficient £ characterises the activity of the actin-myosin network inside the cell, being negative for contractile motors 
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and positive for tensile motors; Sfi is the difference of the chemical potential between Adenosine triphosphate (ATP) 
and hydrolysis products. Let us now fix the boundary conditions: i) the normal component of the velocity at the 
boundary (eq. ([5])) must be equal to the translational velocity of cell along the normal, ii) the other conditions concern 
the mechanical balance of forces on the interface dfl both in the normal and the tangential directions 

<Tnn = -P + 277"^ = CSfi, (12) 

% = v( -7T- + -?r~ ~ \ K \ Vs I = (13) 

\ os on ) 

v n , v s being the velocity normal and parallel to the boundary Sf2 while k is the local curvature. If we introduce a 
tension or capillarity effect, we must modify eq. (|12p and add a contribution proportional to the curvature of the 
interface (Laplace law for the pressure jump). In case of the Laplacian flow with r\ = 0, eq. (flU)) should be ignored. 

For a simple cell geometry as the circular one, the shape equation can be derived by assuming the radial geometry 
for ^ as done in [5]. Nevertheless, sometimes this geometry, which is assumed a priori, cannot allow to satisfy all 
the equations because we face a free-boundary problem, which is especially difficult to analyse for Darcy and Stokes 
flows [6] . Free boundary problems can be solved either by complex analysis [6] or with the help of the Green function 
techniques. The latter case requires numerics while complex analysis is more elegant. Unfortunately, none of these two 
techniques is sufficient to give a solution of the complete problem in the case of eq. (jlip which couples a Laplacian and 
a Bi-Laplacian. One possibility may be that the friction coefficient is large, so that the viscosity effect is located only 
at the boundary layer around the border of our cell. In this case, it will be possible to prove that effects of this viscous 
boundary layer result only in a modification of the boundary conditions by introducing "a capillary term" |12) . 

Let us consider first the case without viscosity: the dimensionless parameter fj — n/(^R 2 ) [5] (R being a radius 
of a circular cell at rest) goes to zero and we can define a velocity potential 3*°) equal to — P/£ . This velocity field 
satisfies the Poisson equation because of the depolymerisation effect. So we have A<P^ — — ka (equivalent to eq. 
and two boundary conditions for the normal velocity (eq. ©) and the other for the normal force balance (eq. (fT2"|) ) 
with r\ = 

v n = v • n = — — = U cos <d - V p , (14) 
an 

P = — (Sfi + capillary terms. (15) 

The capillary terms include e.g. the tension, which exists in the lipid membrane and that the membrane exerts on the 
cell. Note that this problem has some mathematical similarity with the treatment of the sliding drop in the situation 
of partial wetting given by Ben Amar et al. [7] studied with the Schwarz function technique that we present now. 



3 Complex analysis formulation 

Here we develop the complex analysis formulation and Schwarz function technique which is well adapted to bounded 
domains. Assuming the strong friction limit, the Darcy law provides the velocity potential <P^°\ Let us write the 
boundary conditions using complex notations, assuming a pressure jump at the boundary due to capillarity, 

■aT = t - v * (fl> = ia?' (16) 



dn 



n- V<Z> (0) = {Ucosd-V p )+d t n-n, (17) 



where 1? is the angle between the y-axis and the normal n to the interface or the angle between the x-axis and the 
tangent t to the interface. The last term in eq. (flTf takes into account the fact that the border can fluctuate. By 
combining these two boundary conditions, we get 

2 -dr = -d^- l -dT = {-dT- l -^-) e = tV ? e -*2 (1 + e ) + W e ■ (18) 

We define the Schwarz function of the contour ? a such that z = g(z), where z — x + iy is the complex variable and 
z = x — iy denotes the complex conjugate. By using the relation 
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we find that e * = y/g'(z). So finally the complex derivative of the potential along the interface is 

,9(5(0) d f a" \ U 1 

2^ = ^VTM + (^) - i a (i + ✓(*)) + - 2 o t9 - (20) 

Since A^°) = — k d the velocity potential can be written as 

V) = -^2+^(2) + (21) 

where -ff(z) is an arbitrary unknown function of the complex variable z. Combining these two results we find the 
interface equation that we extend to the whole volume of the drop 

- k d g(z) + H'(z) = 2iV P Vfl 7 (M) + i~ (j^j - iU(l + g'{z)) + d t g(z). (22) 

Let us consider the case now where the flow field is along the y axis. As suggested by the experimental results of 
[llj , we decompose the flow into a constant velocity Vo and a flow with non-zero divergence contribution according to 
eq.Q. Then the flow potential becomes 

$M(x, y) = -^y 2 + V y + \{H{z) + H{z)) = *±{z - z) 2 + -|(z - z) + \(H(z) + H(z)), (23) 



which does not modify eq. 

To simplify the notations, we have mentioned only the spatial dependance of the functions <&(°\g(z) and its 
derivatives of first order g'(z) and second order g"(z) but it is obvious that in a non-steady problem they are time 
dependent as well. This equation has been established at the boundary dQ. We extend it to the whole domain fl by 
following the standard strategy for free-boundary problems. Since represents a physical quantity, it cannot have 
singularities inside the domain Q but the Schwarz function and its derivative have singularities inside fi. So we need to 
solve this equation with an arbitrary regular analytical function inside the drop H(z) imposing that g(z) is a Schwarz 
function that is: z = g(g(z)) for arbitrary z inside the drop. We can show straightforwardly that the static disc with 
U = is a solution of the eq. (|2"2")l . Since the Schwarz function of a circle of radius R is g(z) = R 2 /z, g'(z) = —R 2 /z 2 
so y/g'(z) — iR/z, we recover kdR = 2V P . Capillarity gives no contribution in this case. Moreover, it is absolutely not 
obvious to find a traveling wave-solution. Indeed, the leading order singularity related to the velocity U is of the order 
of g' , which has no counterpart in eq. (|22[). and thus cannot be cancelled [7]. It means that such a circular drop cannot 
move with a constant velocity U. Our aim is then to analyse more deeply the polymerisation-depolymerisation process 
of the actin cortex localized below the lipid membrane in order to modify the model which is too much simplistic. Note 
that if the polymerisation process depends explicitly on the polar angle, then the circular moving drop (with velocity 
U along the y-axis) is a solution of our problem. Indeed, by transforming V p , which is the rate of polymerisation, into 
V p (l + acosi?) we get the following flow equation: 

- k d g(z) + H'(z) = 2iV p y/g>(z,t)+ iV p a(l + g'(zj) + - iU(l + g'(z)) + d t g(z), (24) 

and we can identify the drop velocity U as the polymerisation rate times the anisotropy coefficient. Nevertheless, this 
cannot explain the biological processes at the origin of the bifurcation in shape motility. In the next section we analyse 
the stability of such solutions. 



3.1 Linear stability of the static circular solution 



The loss of stability of the static solution may explain the shape bifurcation and the cell motility. The Schwarz 
function technique is very efficient for performing such analysis avoiding complicated geometric calculations. It is why 
we present first the method performed in [5], then the method derived from complex analysis. 

3.1.1 Classical method 

The border is assumed to fluctuate around the circle and its equation is given by 

r = i?(l + e„e f2 "*cosntf). (25) 
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The velocity potential 3>W up to the leading order is 



*(°)(r,0) = --jr 2 + ia„e fi "*r™cosrM?. (26) 



The Laplace equation at the boundary df2 becomes 

?2 



<2>(°)( r ,tf) = -^R 2 - ^ e „ e fi "*cosm? + ^e^iTcosntf = -~(n 2 - l)e n e n - 1 cos m?, (27) 
4 2 4 § it 

which allows to find the relation between a n and e n , given by 

4e„ ^ k d „ 2 7/2 



_ (28) 
Up to the linear order, the normal remains along the radius, so the linear expansion of eq. (I17|) becomes 

- -|-e n e n »* cos nd + Ja„e r2 " t i?™~ 1 cos n§ = Q n Re n e° nt cos nd, (29) 



« B = (n-l)§(l-r 7 7= 5 n(n + l)]. (30) 



and finally the dispersion relation is 

0_ = fr> - 1 1 

i? V V p iR 2 
3.1.2 Stability analysis with the Schwarz function method 

In the presence of small distortions at the contour df2, the Schwarz function is given up to linear order by 

g( z ) = ^(l + e n e^ t (z n + z' n )), (31) 

and H'(z) = b n (n — l)e fi "*z™ _1 . The elimination of the singular mode z~ n in eq. (|2"2"|) gives the following dispersion 
relation 

- k d R 2 e n z' n - x = -V p R(n + IK*""" 1 + ^n{n 2 - I)*-"" 1 + f2 n R 2 e n z~ n - 1 , (32) 

which is exactly the same as in eq. (I30|) . Such analysis shows that the process of polymerisation-depolymerisation 
of actin destabilizes the cell when it is static but the tension in the lipid membrane at the origin of the capillary 
effect may be enough to stabilize the cell for small radius and moderate friction represented by £. Clearly the lipid 
membrane may readjust its tension: it has been experimentally shown that the actin cortex strongly interacts with 
the lipid membrane pQ. 

In summary, this method gives an easy way to prove that a prescribed 2D-shape is a possible solution or cannot be 
a solution to a free-boundary problem. Predictions of more complex shapes than the circle is much more unlikely due to 
the difficulty to satisfy the requirement z — g(g(z)). The algebra for linear stability analysis close to an exact solution 
is straightforward with this method, as shown above. Anyway, we cannot have a motile cell with these boundary 
conditions. Our purpose now is to take into account the anisotropy of the actin network to modify these boundary 
conditions. 

4 Modelling of active gels as two-dimensional ordered fluids 
4.1 How to introduce the anisotropic nematic tensor 

We describe the actin filaments as nematic objects with a local averaged orientation represented by a unit vector p 
called the director varying from point to point inside the cell. The order inside the cell is measured by the orientational 
order parameter q. When q — the filaments are randomly oriented and when q = 1 they are perfectly ordered. 
Therefore, similar to nematic liquid crystals, actin filaments may be described by a tensorial order parameter [131114] 

Q = ?fp®P-^J or Qij = q(PiPj ~ <W 2 )) ( 33 ) 
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where I is the two-dimensional identity tensor. One should remember that Q is a two-dimensional traceless symmetric 
tensor. Moreover, there is an intrinsic symmetry of the problem, namely the tensorial order parameter Q, associated 
with negative q and director p, coincides with the one associated with positive q and director orthogonal to p. Thus, 
without loss of generality, we restrict our analysis to non-negative values of q. 

The hydrodynamic equations for the order parameter Qij of actin filaments and velocity of the gel can be regarded 
as an extension of the continuum theory for nematic liquid crystals with tensorial order [151116] towards the theory of 
active gels i8'..9lH0]. In the cell cytoskeleton, the molecular motors (mostly myosin II) move along the actin filaments and 
modify the organization of the actin network [T] allowing the deformation of the cell and its possible displacement on 
a substrate. These active motors transform the chemical energy coming from Adenosine triphosphate (ATP) hydrolys 
into mechanical energy generating a nonzero active stress af?* = — C<fyt(Qij + <%)• According to the sign of £, negative 
(positive) the stress acting on the actin network is contractile (tensile). The coupled equations of motion for the order 
parameter Q and the velocity field v are given in a macroscopic approach by 15,8IIWIT)] 

£v i =a i j ! j, (34) 

—f^T = Qij,kVk = PiUij + \-Hij - XSfiQij, (35) 
jjt pi 

where represents the stress tensor, Vij the velocity gradient and Uij the deviatoric part of the velocity gradient: 
uij = (vi t j + Vj j i)/2 — (vk,k/2)$ij- The operator D/Dt denotes the convective time derivative. Furthermore, the 
summation convention has been assumed. 

The stress tensor consists of four contributions, an = o\ s „ + cr" + ap/ + of/* where 

<, l'<\, I -*'/'-..,. •'://,,. <J a = H ik Q kj :-Q ik H kj , (36) 

ff tf r = "Jt-Gmj, <r act = -CWtf • (37) 

where P is the hydrostatic pressure, r\ the ordinary viscosity mentioned before, j3% the rotational viscosity and W the 
free energy density. The simplest form of this free energy density includes the elastic nematic energy within the the 
one-constant approximation and the Landau-de Gennes potential: 

W = A'|VQ| 2 + aTr(Q) 2 +cTr(Q) 4 , (38) 

where K and c are positive constants, whereas a changes sign at the temperature of the nematic-isotropic transition. 
The molecular tensor field H is the opposite of the functional derivative of the free energy density W with respect to 
the order parameter Q 

SW ( dW \ dW 



H v = -JrT= -*rT- ( 39 ) 

In the following we neglect afj and tr^ r , because active gels are usually supposed to be poorly ordered, namely |Q| <C 1, 
and thus only the lowest order terms are taken into account in eq. (|34p . 
The tensorial equation f|35[) can be splitted into two scalar equations 



V • V? = -j- ( K(Aq - 4q\Vp\ 2 ) -aq-^q 3 )+ 2/? lP • (Vv)p - ft divv - \6fiq, (40) 



p 2 V v ' 2 

qp X (Vp,v) = p x 



f [gAp + 2(Vp)(Vg)] +ft(Vv)p 

P2 



(41) 

The presence of active motors, represented by A5/i, leads to the renormalisation of the term linear in q. 
4.2 Circular geometry 

Let us first consider these equations in the polar coordinate system {r, tp}. The Laplacian of the scalar order parameter 
q becomes 

Aq = -d r {rd r q) + ^d w q. (42) 
We parametrise the director as p = cos 8e r + sin 8e v and so we have 

Vp = p rr e r (g) e r + p rv e r ® e v + p vr e v ® e r + p vv e v ® e v , (43) 



8 



M. Ben Amar et al.: Cell motility: a viscous fingering analysis of active gels 



where 



p rr = d r (cos 9), 
p vr = <9 r (sin0), 



|Vp| 



Pr 



Prtp ' 

Ptpxp 
„2 



(dip cos 9 — sin t 



IPr^p Pipr Pifip 



- (d v sin ( 

r 

(d r 9) 2 4 



COS ( 



We consider first a velocity field in the general form v = v r e r + v v e v , yielding the velocity gradient tensor as 



1 



and thus 



Vv = d r v r e r ® e r + — I d r v v H — (d v v. 



1 



(e r 



1 



e r ) + -(d<pVu, + v r ) e ¥ 



divv = d r v r + - {d v v v +v r ), 



p • ( Vv)p = cos 9d r v r + sin 9 cos 9 ( d r v v H — (d v v r — v v ) ) + 



sin 2 9 



{d v V v + V r ), 



p x (Vv)p 



cos 29 ( d r v v H — (d v v r — v v ) \ + sin 29 ( —d r v r H — {d v v v + v r ) 



(44) 
(45) 
(46) 

3 

(47) 

(48) 
(49) 
(50) 



Then eqs. gOj) and (gTJ take the form 

KAq - AqK ( (d r 0) 
+ Pl 



v r d r q + -v^d^q = \- 
r p 2 



^(9 V 9 + 1) 2 



q(a + \8nP2) 



2 q 



+ 



1 



cos 20 d r v r (d<pV v + v r ) I + sin 26* d r v v H — {d v v r — v v ) 



(51) 



and 



q(v r d r 9+-^(d v 9 + l)) = 



KqA9 + 2K ( d r 9d r q 
+ 4 



rd v q(d v 9 + i; 



1 



cos 20 d r v v + -(d^pVr —v tp )\+ sin 20 —d r v r + -{d v v v + v r ) 



1 



(52) 



Both equations are essentially non-linear and the orientation of actin filaments is coupled to the flow, which will 
be found at the last section [6] With the flow equation, we should deal with three coupled non-linear partial differential 
equations (P.D.E) for q, 9 and 'f', the stream function. There is no hope that these equations can be solved without 
assuming a certain relationships between the coefficients as well as the form of flow. At the boundary, the cell membrane 
together with immersed proteins imposes a preferred orientation on the actin filaments. Therefore, we will assign at the 
boundary a fixed angle 9q between the normal and the fiber director p and a fixed degree of orientation q, which may 
depend locally on the shape of the membrane via its curvature. In a narrow region near the boundary we suppose that 
q and 9 change rapidly. Then the solution can be written in power series of e, which is a small dimensionless parameter, 
relating the width of a boundary layer, defined below, to the typical size of the cell R. To formulate the boundary 
layer problem, one needs to take into account non-linearities, and one expects to find the soliton-like behaviour of q 
and 9. Far away from the boundary, at the core of the cell, we have two cases: either isotropic state with 'melted core' 
(qo = 0) or the uniformly aligned core with degree of orientation given by 



2(a + X6^ 2 



> 0. 



(53) 



These two cases are shown schematically in the fig. [5J Depending on the activity of motors XSfj, (positive or negative) 
the sign of a + \8fiP2 can change. We will treat separately the cases with a melted core and an aligned core. 



5 The bulk and the actin cortex 

5.1 Melted core: boundary layer for q 

Let us first consider a simple case, assuming that the angle 9 is constant in the small region near the border, and 
its value is imposed by the cell membrane, so that 9 — 9q. Given that actin filaments approach the lipid membrane 
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A) Melted core 



B) Aligned core 




Fig. 2. Depending on the activity of motors A<5/x we can have two cases with A) melted core (a + A5/i/?2 > 0) or B) aligned 
core (a + XS^fa < 0). In the case A) we define the width of a boundary layer (ei-R) for the change of the degree of orientation 
q; in the case B) we define the boundary layer (eR) for the change of the director p. 

almost symmetrically, with some preferred angle either ~ 20° or ~ —20° [U[S]j the average orientation at the border 
is along the normal, thus 9q — 0° and the nematic tensor Q is 



Qrr Q<p<p ^ ' and Qrtp Q tpr 0* 



(54) 



Furthermore, we assume that the velocity field v — v r e r + v v e v is small and can be disregarded for a moment (below 
we define a small parameter t\ and the order of smallness for v r and v<p). Thus, eq. (|51[) for the degree of orientation 
takes the form 

(55) 



KAq - q{a + XSfxfo) - -q 3 = 0. 



Based on this equation we can introduce the dimensionless parameter 

1 / 



K 



ei 



(56) 



i?Y i a + Afy&r 

relating the thickness of the boundary layer to the typical size of the cell R. In the case of liquid crystals, t\ is the 
ratio of the nematic coherence length y/ K/a, which is of the order of nanometers, to the characteristic length of the 
system R ~ 1 — 10 /xm, yielding ei < 10 -3 . We can now estimate the strength of the advection term v r d r q or v r d r 9. 
Assuming that the main flow is given by —kdR/2, the dimensionless parameter ^kd/a is of order e\ . In a thin region 
near the boundary, we suppose that q can be written in power series of a small parameter ei 



q = q {0) + e ig W 



Let us introduce the rescaled variable 



i? 



= 1 - 



(57) 
(58) 



measuring the distance from the boundary. So ? = corresponds to the border of the cell r = R, and <; — > oo 
corresponds to the interior border of the boundary layer (see fig.[2K)- At the leading order eq. (|55|) takes the form 



<W 0) ~ q {0) - w(<7 (0) ) 3 = 0, where u 



2(a + \Sfip 2 ] 



>0, 



(59) 



which should satisfy the asymptotic conditions 



7 (0)| 



3c<7 (0) l 



(60) 



and the boundary condition at the border of the cell 



(61) 
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Multiplying both sides of eq. ([55)1 by d^q^ and integrating we find 

(9fg (0) ) 2_ (g (0) ) 2_| (g (0) ) 4 = 0! (g2) 

where the integration constant was chosen so to satisfy the asymptotic conditions The solution of this equation 
admits the form 

g (0) (0= ■ f^ ry C = a S inh(^M. (63) 
smh^ + G) \ q J 

As expected, in the leading order we found a soliton-like behaviour for q. In the case of a non-circular cell, a local 
boundary layer analysis must be done with a coordinate frame defined from the normal and the tangent at the point 
of the interface Oil considered. Then s is simply the distance taken along the normal modified by the local curvature. 
The set of coordinates is then <; and the arclength s. It turns out that for the leading order the equations are quite 
the same as the one derived from the cylindrical coordinate frame with some obvious adaptation which will be given 
without demonstration |12j . In a local curvilinear coordinates the Laplacian up to the next order is given by 9 SS — kc^, 
where k is the local curvature, k = dd/ds = —dip/ds, chosen to be negative for convex interfaces. The components of 
velocity field v rj v v being of the order of e 2 , the next order term in the expansion eq. (1571) should satisfy the following 
equation 

d^qW -qM^l + S^f) = + (64) 

For convergence of the series, expansion represented by eq. (|57|) . we look for a solution of eq. (|64)) with vanishing q^ 
for ^ = and q — > oo. The adjoint theorem gives us the value of d^q^ at cr = for such a solution, an information 
which is useful for future development of the analysis. 



5.2 Aligned core: boundary layer for p 



Now, the core of the actin network is ordered (see fig. HJ3). The coupling of the director's orientation with the flow, 
given by the term proportional to /?i in eqs. (jSTj) . (l52j) . leads to the alignment of actin filaments in the direction of the 
flow. Without loss of generality we assume the flow to be in y-direction, as in section[5] Therefore, the core of the cell is 
expected to be aligned in the same direction as shown in fig.[U The main velocity field is of the form v = (Vq — kdy)e y 
with divv = — kd (see section [5] and eq. ([3])). Note that at this stage, we do not have information on the value of Vq. 
Then eqs. (J5TJ) and ([52]) take the form 



(Vq — kdX sin <p) ( sin Lpd r q H — cos ipd^q 



1 

~2 



and 



KAq - 4qK ((d r 6) 2 + ^{d v 9 + l) 2 ^j -q{a + \5[i(3 2 ) - °-q z 



+ /3ifc d cos2(^ + 9), (65) 



q(Va — kdX sin ip) ( sin (pd r 9 H — cos <p(d v 9 + 1) 

1 



KqA9 + 2K[d r 9d r q + —d v q{d v 9 + 1 



-Ay sin 2 (<p- 



9). (66) 



The boundary layer thickness (ei ) for the degree of orientation was defined in eq. (|56l) , and it was a function of the 
activity of molecular motors XSfj,. We may have different situations with a double boundary layer structure, depending 
on the values of the coefficients in eqs. (|d3|) and (pp| . However, it is reasonable to assume that near the boundary the 
scalar order parameter q changes more rapidly than the director. Thus we assume in this section that q is constant, 
given by the equilibrium value in the bulk q = q (eq. (j53|) ) and qo/(3i is of the order of e 2 , the dimensionless thickness 
of the boundary layer for the angle 9, defined as follows 



qoK 
\Mk d ' 



(67) 



The transport coefficients /?i and 02 can be written in terms of the rotational viscosity 71 > and the torsion coefficient 
72 as Pi = —9072/71, P2 = 7i/(2g 2 ) [TMT4"] . The torsion coefficient 72 = r\\—r)i [H], where 771 is the Miesowicz viscosity 
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when p is parallel to v and 772 is the Miesowicz viscosity when p is parallel to Vv, can be positive (e.g. HBAB in 
nematic phase or negative (e.g. MBBA in nematic phase [12]). Since we do not know neither the sign nor the 
value of 72 for actin filaments in cells, we used the absolute value of pifii = —72/(2(70) in eq. (1571) and in the following 
assume 72 > 0. Taking into account that e\ < e (eqs. (1551) . (|5T|) ). the following relationship between the coefficients 
holds 

ei<e ^> \a + X6 t i(3 2 \> \(3i(3 2 \k d /q a . (68) 

To estimate the value of e we take the elastic constant for liquid crystals K ~ 6-10 -12 N and assume that |/3i/32| — 72 — 
77, which is a poorly measured value, varying between ij ~ 0.1 Pa-s |18j . r] ~ 300 Pa-s [19] and r\ ~ 10 3 — 10 5 Pa-s |10) . 
Assuming the r\ ~ 10 4 Pa-s and kd ~ 0.2 1/s [5] we find e ~ 0.01 for the cell with radius R = 10 /im and e ~ 0.1 for 
R = 1 /xm. The smallness of e justifies the assumed boundary layer approximation, however, the lack of information 
about the value of the torsion coefficient 72 and the rotational viscosity 71 remains a sensitive issue for our model. 
We aim at solving the boundary layer problem for the director, by expanding the angle 9 as a perturbative series 

= 0<°> +e# (1) + e 2 (2) + ... (69) 

Let us introduce the rescaled variable p 



r 
R 



1 - ep, (70) 



which measures the distance from the boundary, so that p = corresponds to r = R and p — > 00 corresponds to the 
interior border of the boundary layer (see fig. HJ3). We should satisfy the boundary conditions for the orientation at 
the zero-order approximation, so that 



(o) Uoo = J-^ #U=«o, (71) 

whereas for the higher orders 

9 (i) U = « (i) lH0 = <l. i = 2,3,... (72) 
The leading order approximation for the orientation satisfies the following equation 

2cV6»(°> + sm2(tp + 6^) = 0, (73) 

with the first integral given by 

2^/)+ S m2(H^) = 0, -> (d^r- 0082 ^ 9 ^ =C{ V ). (74) 

The constant of integration C(ip) can be found from the condition in the bulk eq. (jTT]) . yielding C(ip) — 1/2 and 
consequently 

(d p 9^) 2 = cos 2 (^ + (o) ). (75) 



The further integration gives 

n (a\ „ M(^)e p -1\ \ 1 + tan^ 
6» (0) = + 2 arctan — ^ , A(<p) = s- 2 — , (76 



where the constant is chosen to satisfy the boundary condition at the membrane p °\ = 0q (eq. (jTTj) '). 

To illustrate the convergence of the series we need to find the next order term in the expansion. By substituting 
eq. (|69p into eq. (|66[) . taking into account the contribution from the LHS, and collecting the terms proportional to e, 
the equation for 8^ takes the form 

d pp 6V + 90.) cos 2fo + = -«0„*<°> f 1 - + . (77) 

V \Pl\ k d \pi\R ) 

V v ' 

=V(v) 

Substituting the results for 9^°', eq. (f?6"|). in eq. ([77]) we find the equation for 



2Ae p 



1 + ^e 2 " 



2 



- 1 



1 + A 2 e 2 » 



Vfo), (78) 



12 



M. Ben Amar et al.: Cell motility: a viscous fingering analysis of active gels 




0(0) 

0(1) 



10 



Fig. 3. The plot of 6> (0) and 8 m given by eqs. (JT6J) and (j 79f) as function of p for ip = 6q = and k = —1. 



and the solution in a compact form 



(1) AVp + sinhp 



(79) 



which satisfies the boundary conditions flW (0) = (oo) = 0. In fig. [3]we compare the first two terms of the expansion 
of 9 (j6"9")l and notice that the amplitude of 9^ is smaller than in the whole range of p. Nevertheless, eq. (|T9")) is 
valid if and only if q is kept constant up to order e 2 in this boundary layer. This restriction which is absolutely not 
necessary deprives us of a degree of freedom which would be very useful for the next orders of the asymptotics when 
the flow field and the shape of the cell is considered. A much better strategy will be to solve eq. and eq. (pjf by 
keeping the expansion eq. (1691) for 9 but the following expansion for q 



q = q {l + eq^). 



(80) 



This viewpoint does not modify our leading order and we consider now the modification of the hydrodynamic equations 
induced by the actin cortex. 



6 Flow in the actin cortex 

Actin cortex exerts forces on a lipid membrane, resulting in the modification of the boundary conditions (eqs. (|16p and 
PT| ) and consequently the shape of the drop, eq. (|2"2")l . Let us consider the balance of forces given by eq. Since 
we are in the limit of the boundary layer we adopt a Cartesian system of coordinates. Later on, we will change to a 
local coordinate system more suitable to fix the boundary conditions. Taking the divergence of this tensorial equation, 
we find a modified Stokes equation (see eq. ©), 

BP BA- 

where 

Aij = -(SfiQij - 2Px(KQij^k - aQij - 2cQijQ k iQki), (82) 

is the stress tensor (see eqs. (155]) . (|57| and (|59")l ), accounting only for the lowest order contributions as was referred 
above. Introducing the stream function \P as before in eq. (O, so that 



we can rewrite the eq. (181 1) in the form 



m BP 

*7* + 7k = " A 7* +VA l" (84) 

d&\ BP B$ 
-e(A W + % + _j + _ = _^_ + V A| 1 , l (85) 
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A) k d = V = 



B) k d + 




C) Vo ^ 



Fig. 4. Flow in the actin cortex for (v x , v y ), given by eq. (|83jl for the stream function IP (eq. !94p only (A) and for IP supplemented 
with a non-zero value of kd (B) or a non-zero value of Vo (C), which are the dominant contributions. The meaningful part of 
the flow is only near the border, where the actin cortex is localised. 




where VA|^ = Ayj. Taking the cross derivatives of the above equations (<9 y eq. ([51)1 — c^eq. (|85[) ) we eliminate the 



pressure P and find 



£AV = r,A 2 & + d y VA\ x - d x VA\ t 



(86) 



This equation will replace the hydrodynamic equation One can also recognize the part of eq. (fTTj) derived in the 
section [2l As previously, we neglect the viscosity as was already stated in section [2l so i] — 0. Note that in this limit 
far away from the actin cortex Alp^oo ~ 0, one recovers the Laplace equation for the stream function and therefore 
the velocity potential formulation considered in the section [2] 

Coming back to the polar system of coordinates, so that {x, y} — > {r, ip\, we have 



VA| r 



dA rr 1 / dA rif 

dr 
OA 



r \ dtp 



ipr 



dr 



OA 



rr 



yielding eq. ((86]) in the case r/ = as 
ld 2 A rr i S* 2 



\_(o- \ 



dA r 



dA^ 



dtp 



ld_ 

r dr 



A 

-A 



A 



'r'r 



A 



pA vr 

dr 



ld 2 A u 



drdip 



r drdtp r 2 \ dtp 2 dtp dtp 
Substituting the equilibrium value of the scalar order parameter go (eq. (1531) ) into eq. 

A i:j = -((Sfi + 2/3 1 /3 2 A^)Q„- - 2^KQ ljM . 



1 dA,p r 
r dr 



I dA r 



dr 



we get 



(87) 
(88) 

(89) 
(90) 



This implies that the Laplace law of capillarity is modified by adding A rr . Next, based on eq. (|68p we compare two 
terms in eq. ((9"01) , and find that when X8fi ^> kd we can account only for the active stress and omit the last term which 
is elastic. Indeed, according to experimental observations [TT| . the motility of cell can be related to the increase of the 
activity of the motors, therefore it is reasonable to assume that XSfi increases, yielding only the term linear in Qij 
within our boundary layer approximation 



(91) 



We will solve eq. (|89[) for the stream function within the boundary layer approximation. Since in section 0] we 
assumed the leading order contribution to the flow in the form v = (Vo — hay) e y , the stream function <?", representing 
here the higher order corrections, takes the form 



(92) 



Substituting this expansion together with eq. (|9Tj) into eq. ((89)) and introducing the parameter £ = ^/(C5/x+2/3i/?2AJ/i), 
so that £ ~ 0(l/e 3 ), which is indeed the limit of strong friction assumed above, we find 



£(d PP ^ (1) + ed pp *W 



e K d p & w 



d n(o) +£( 9 n(i) 

u pp^rtp ~ cu pp^rip 



260 ^psQrr > 



(93) 



14 



M. Ben Amar et al.: Cell motility: a viscous fingering analysis of active gels 



where we used the fact that Q is a symmetric and traceless tensor, and we have introduced a local curvilinear system 
of coordinates as above with the local curvature K = dd/ds = —d(p/ds, negative for convex interfaces. The leading 
order contribution to the stream function is 

= ^ = ^sin(2^°)), (94) 

where qo and 9^ are given by eqs. (|53p and (|76p . respectively. Taking into account this correction to the main velocity 
field v = (Vq — kdy)e y , considered in the section [H in fig. [4] we visualise the flow given by eq. (l83|) . The single 
contribution of the stream function IP ((eg. fig. UK) suggests that the flow near the border, where the actin cortex 
is localised, has a preferred y-direction. Therefore, we may expect to find the motility of the cell as will be shown in 
the next section. In the case of the melted core Q rip = (eq. ([54]) ) and we find 

G^ (2) 2 dQ rr 

~ 5 = Q ' 95 ) 

dp £ ds 



7 Boundary conditions 

Due to the Cauchy relations 

n • W = -t • V<F, t • \7<P = n • V<f , (96) 
we can obtain the velocity potential from the stream function 

$ = <|>(°) + e 2 $^ + + eV 3 > + . . . , (97) 

where is given by eq. (|2U)) used in the section [2] to formulate the boundary conditions eq. (|2"Tj) and being of order 
e 2 . The next order contribution to the boundary conditions is calculated by taking into account eq. (|94p. resulting in 



0$(i) _ ggW _ | K | g0 9sin(2gW) 



dr 9s 2£ dip 

d$W g&W q <9sin(26>(°)) 



= 0, (98) 

p-s-0 



ds dr 2R£ dp 



= --^cos((9 + V5)cos(2(9o). (99) 



To sum up we have obtained a modification of the Neumann boundary condition: <l>^ = —q/(R£) which has some 
consequence for the melted core solution if and only if the prescribed value of q at the boundary depends explicitly on 
the arclength. In the absence of biological information, we have taken this value constant along the cell boundary and 
this contribution does not modify the free boundary problem. In the case of the oriented core we have a modification 
of the Neumann condition and we can evaluate its consequence on eq. (|2"2"|) . Cancelling the value of #o (averaged value) 
and noticing that <p has to be identified to 7r/2 — § the left- hand-side of eq. (|I6I) has the following additive contribution 
given by 



£ - T^b -7= - vVW • (100) 



As a consequence the left-hand side of eq. (|22[) is also modified by a term which is — iqo{l — g'(z))/(R£) which 
compensates exactly the singularity induced by the velocity and we derive: 

U=-^ = -^CSp + 2p 1 p 2 XSp)j. (101) 

First we note that the velocity U is positive for contractile motors in agreement with experiments. We can roughly 
estimate the value of this translational displacement of the cell. Taking the averaged orientation of actin filaments at 
the border as zero, 8q ~ 0, the radius of the cell R ~ 1 — 10 p,m, the value for the measured friction coefficient in 
keratocytes £ ~ 10 5 Pa-s//z 2 [5J, and the activity of motors (5p ~ — 10 3 Pa, corresponding to contractile motors [5 
we get an estimate of the first term in eq. (|101j) as U ~ 1 — 10 nm/s. We also assumed a perfect ordering in the center 
of a cell, go = 1 in eq. (1531) . The second term can be rewritten taking into account that 2/3i/?29o = — 72, which gives 
\8 / {R£) ■ The viscosity of actin filaments r\ ~ 10 4 Pa-s [5lll0j. which is orders of magnitude higher than for any 
nematic liquid crystal. Assuming 72 ~ r\ and XSp ~ 5kd ~ 1 s -1 , we find the velocity U ~ 10 nm/s of the same order 
as above. 
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8 Conclusions 

We have shown that it is possible to explain the cell motility through the reorganisation of the actin network in a two- 
dimensional cell in the limit of strong friction interaction with a substrate. The cytoskeleton dynamics is described 
by the hydrodynamics equations for active gels. Due to polymerisation occurring below the lipid membrane and 
depolymerisation at the other end of the filament, a flow of F-actin filaments exists that we modelled by a Darcy 
flow. The fact that the organisation of the cortex imposed by the membrane is limited to a short distance from the 
border allows to treat the equations of active gels for an arbitrary cell geometry. Nevertheless, to leading order in the 
thickness boundary layer, we have shown that a circular cell with no order orientation for the F-actin filaments in the 
bulk is static while an averaged longitudinal order induces an hydrodynamic flow compatible with the model (bulk 
hydrodynamic and boundary conditions). Modification of the cell shape can be obtained by extending the present 
work to next orders and by introducing biological informations necessary to improve the description. 
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and Jean-Frangois Joanny. O. M. was supported by the French National Research Agency (ANR), grant ANR-07-BLAN-0158. 
G. N. was partly supported by University "Pierre et Marie Curie". 
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